#!/usr/bin/perl -w
use strict;
package interval_set_mark_overlap;

require Exporter;
use strict;
use warnings;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(	
					interval_set_mark_overlap
				);

sub min($$)
{
	return ($_[0] < $_[1]) ? $_[0] : $_[1];
}
sub max($$)
{
	return ($_[0] > $_[1]) ? $_[0] : $_[1];
}

# takes two sets of intervals (begin and end, begin < end)
# and marks those that overlap
sub interval_set_mark_overlap(\@\@\%\%)
{
	my ($seta, $setb, $overlapa, $overlapb) = @_;
	return if (!@$seta || !@$setb);
	my ($posa, $posb) = (0,0);
	while (1)
	{
		# if a is before b, increment a until end of a
		if ($seta->[$posa][1] < $setb->[$posb][0])
		{
			last if (++$posa >= @$seta);
			next;
		}

		# if b is before a, increment b until end of b
		if ($setb->[$posb][1] < $seta->[$posa][0])
		{
			last if (++$posb >= @$setb);
			next;
		}

		# overlap
		{
			my $max_start = max($seta->[$posa][0], $setb->[$posb][0]);
			my $min_end = min($seta->[$posa][1], $setb->[$posb][1]);
			my $overlap_len = $min_end - $max_start + 1;
			$overlapa->{$posa}[0]+=$overlap_len;
			$overlapb->{$posb}[0]+=$overlap_len;
		}

		if ($seta->[$posa][1] < $setb->[$posb][1])
		{
			last if (++$posa >= @$seta);
			next;
		}
		else
		{
			last if (++$posb >= @$setb);
			next;
		}
	}

	#
	# 	calculate percentage coverage
	# 
	for my $pos(keys %$overlapa)
	{
		my $overlap = $overlapa->{$pos}[0];
		my $len = $seta->[$pos][1] - $seta->[$pos][0] + 1;
		$overlapa->{$pos}[1] = $overlap / $len;
	}
	for my $pos(keys %$overlapb)
	{
		my $overlap = $overlapb->{$pos}[0];
		my $len = $setb->[$pos][1] - $setb->[$pos][0] + 1;
		$overlapb->{$pos}[1] = $overlap / $len;

	}
}

sub test()
{
	my @a = ([1,300], [501, 700], [ 801, 900 ], [1101, 1300]);
	my @b = ([121,220], [451, 800], [1300, 1499]);
	my (%a_results, %b_results);
	interval_set_mark_overlap(@a, @b, %a_results, %b_results);
	print STDERR "[".$_->[0], ",", $_->[1]."]," for (@a);
	print STDERR "\n";
	print STDERR "[".$_->[0], ",", $_->[1]."]," for (@b);
	print STDERR "\n";
	
	print STDERR "[".$a[$_][0], ",", $a[$_][1]."]=$a_results{$_}[0]:$a_results{$_}[1]," for (keys %a_results);
	print STDERR "\n";
	print STDERR "[".$b[$_][0], ",", $b[$_][1]."]=$b_results{$_}[0]:$b_results{$_}[1]," for (keys %b_results);
	print STDERR "\n";
}
#test();
